陈 瑜，张小艳，许汉华，刘文连，李 泽

(1.昆明理工大学 建筑工程学院，云南 昆明 650500；2.昆明理工大学 电力工程学院，云南 昆明 650500；3.中国有色金属工业昆明勘察设计研究院有限公司，云南 昆明 650051)

作为世界上地质灾害最严重的国家之一，我国的地质灾害不仅种类繁多、难以监控，而且分布范围广、破坏性强。滑坡是我国地质灾害的主要类型，特别是近年来西南地区滑坡灾害逐渐严重，对人民群众的生命财产安全造成了严重危害。地震作为岩质边坡失稳破坏的主要影响因素，其作用下的岩质边坡动力稳定性是保证工程安全建设和运行的关键。要准确计算岩质边坡的动力稳定性需同时考虑地震波传播的时间效应和空间效应。近年来，塑性极限分析法被广泛运用于边坡稳定性的分析之中，其中众多学者运用塑性极限分析法对边坡静力、动力稳定性分析进行了卓有成效的研究工作[1-5]。但许多关于边坡动力稳定性分析问题尚有待于进一步研究。

目前，岩质边坡动力稳定性分析方法主要包括：拟静力法[6-7]、拟动力法[8-10]、基于振型叠加的反应谱法[11-12]以及动力有限元强度折减法等[13]。这几种方法均有其优点，但也有其不足之处。拟静力法作为一种简单高效的地震模拟方法，被广泛运用于各个工程领域中[7]，拟静力法可考虑边坡地震拟静力荷载随边坡高度的变化，但其不能考虑地震的“行波效应”且放大了地震荷载效应；拟动力法作为拟静力法的一种改进，可同时考虑地震波在岩体中传播的时间效应，但由于其将地震波假设为简谐波导致计算结果与实际地震所产生的结果有一定差异；反应谱法基于动力频域分析，虽然在结合有限元方法等方面有优势，但不适用于对时域分析。另外，拟静力法和反应谱法均只能得到单一的边坡安全系数值，而不能获得边坡安全系数的时程曲线。

在当前地震作用下的边坡稳定性研究分析中，大多数仅以某一时刻的安全系数来评价整个边坡的稳定性，这种分析方法显然不够全面。为了进一步完善对岩质边坡的动力稳定性的研究，应增加动力稳定性的时程分析方法，以弥补现有方法的不足，并为地层岩质边坡的动力稳定性研究提供依据。时程分析法作为一种高效的动力分析方法，不仅综合考虑了地震动各个要素对研究的影响，而且能够分析研究对象的能量损失情况，是近年来许多国家抗震设计规范或规程的分析方法之一[14]。目前时程分析法已被推广运用于诸多工程领域，如地下深埋结构相互作用的时程分析[15]、海底管道抗冲击性能的时程分析[16]、T型换乘车站结构抗震设计的时程分析[17]、地下结构抗震性能的时程分析[18-19]、抗爆结构的动力响应分析[20]等。

为了解决岩质边坡稳定性分析中存在的以上不足，本文将塑性极限分析法与时程分析理论相结合，提出岩质边坡稳定性的时程分析下限法。

本文将岩质边坡岩体当作以岩块单元和结构面组成的系统，本系统带有非连续介质力学的典型特征。系统中较为完整的岩块强度比较高，可把滑面上方的岩体当作岩块进行离散。本文为简化计算进行如下假设：(1)滑面上方岩体的破坏在相邻岩块结构面上发生；(2)岩体之间的结构面可以当作理想状态下的刚塑性材料；(3)结构面存在发生剪切破坏的可能性；(4)岩块与岩块间存在相对滑动的可能性。

为了更好地体现地震作用下岩体的移动情况，本文采用岩块单元来离散岩质边坡，总共可将岩质边坡离散为8个块体单元，滑动面以上部分被离散为5个块体单元。岩块单元和相邻岩块结构面的受力情况如图1、2所示，图中(x,y)表示总体坐标系，(nk,sk)表示局部坐标系。

图1 地层岩质边坡的岩块单元i的受力分析图

结合图1、2，岩块所受的各作用力如下：

岩块i形心gi处的体积力Fia：

Fia=(0,Gi)T(i=1,2,…,Ns)

(1)

式中：Ns为滑动部分离散所得岩块的总数；Gi为作用于岩块i形心gi处的自重力，N。

作用于岩块i形心gi处的等效外荷载力Fib：

Fib=(fxi,fyi)T

(2)

式中：fxi为x方向的等效外力，N；fyi为y方向的等效外力，N。

作用于相邻某两个岩块间结构面形心处的内力向量Qk：

Qk=(Nk,Vk)T(k=1,2,…,ne)

(3)

式中：Nk为作用于结构面k的形心Pk处沿外法线nk方向的法向力，N；Vk为作用于结构面k的形心Pk处沿切线sk方向的切向力，N；ne为整个滑动部分离散后结构面的总数。

作用于岩块i形心gi处的地震惯性力向量Fie：

Fie=(Fhie,Fvie)T

(4)

地震是岩质边坡发生破坏坍塌的主要原因之一，本文以拟动力法相关理论为基础模拟岩块受到的地震力，岩块的地震加速度包括水平和竖直两个方向，可分别按照下述方法进行计算。

(1)岩质边坡块体的水平方向地震时程加速度按下式计算：

(5)

(i=1,2,…,Ns)

(2)岩质边坡块体的竖直方向地震时程加速度按下式计算：

(6)

(i=1,2,…,Ns)

岩质边坡滑动面以上的部分被离散为岩块单元和结构面的系统后，以极限分析法理论为基础建立边坡的静力许可应力场，岩块单元必须同时符合对应的平衡方程与屈服条件。本文进行极限分析时以岩质边坡稳定性的安全系数作为目标函数，以结构面的内力作为决策变量，在同时满足平衡方程和屈服条件等约束条件的情况下，建立岩质边坡动力稳定性的非线性数学规划模型。

把岩质边坡的安全系数作为目标函数，进而求解边坡安全系数的最大值，目标函数如下：

maximize:Kj(j=1,2,…,Nt)

(7)

式中：Kj为第j个时步岩质边坡动力稳定性的安全系数;Nt为地震振动时长内的时步数量；maximize为“使之最大”。

若图1所示的岩块单元在受到各种力作用的情况下仍可保持平衡，则考虑容重超载的影响，岩块单元i的平衡方程如下：

(8)

(9)

如图2所示的结构面k在局部坐标(nk,sk)中存在有作用于形心上的内力向量Qk，因此结构面上存在剪切破坏Mohr-Coulomb的屈服条件：

图2 地层岩质边坡相邻岩块单元i与j的受力分析图

(10)

式中：Nk为结构面上的法向力，N；Vk为结构面上的切向力，N；lk为结构面k的长度，m；ck为对应结构面k上的凝聚力，N/m；φk为对应结构面k上的内摩擦角，(°)。

根据以上条件可以得出考虑岩块平动时结构面的广义屈服条件，如下所示：

CkQk≤Bk

(11)

(12)

(13)

以极限分析的下限定理为标准，可以推导出岩质边坡的静力边界条件：

(14)

将目标函数式(7)、块体单元平衡方程式(8)、结构面屈服条件式(11)以及静力边界条件式(14)相结合，可以得到地震作用下岩质边坡动力稳定性的线性数学规划模型，如方程组(15)所示。

(15)

式中：mi为第i个块体的质量，kg;Ts为地震震动的总时长，s；Nt为地震震动时长内的时程数量。

本文经过一系列的分析求解后建立了地震作用下岩质边坡动力稳定性时程分析的极限分析法数学规划模型，本模型是典型的非线性规划模型。针对所建立的非线性规划模型，使用Matlab编写了对应的计算程序。对本程序优化后，求解可获得岩质边坡安全系数的下限解，求解流程如图3所示。

图3 岩质边坡安全系数非线性规划模型的求解流程图

本文算例1选取一个人工开挖后的岩质边坡，如图4所示。表1为此算例边坡的岩体材料物理力学参数表。算例1边坡高度为30 m、坡顶宽度30 m、坡角为67.0°；边坡体内滑动面H1、H2的倾角分别为35.0°、37.5°；岩体节理L1、L2相互平行，两者的间距为14.9 m，倾角均为83.0°。岩体被划分成8个岩块。本文将材料的容重设为定值，仅考虑结构面凝聚力和内摩擦角的变异性。岩质边坡内无地下水作用。

图4 算例1岩质边坡的几何形状示意图(单位：m)

表1 算例1岩体材料的物理力学参数

本算例采用1988年11月澜沧耿马地震实测得到的地震波，地震加速度峰值按7度地震选取，为0.2g。在边坡坡角处输入的水平向、竖直向地震加速度时程曲线如图5所示。地震波间隔0.02 s输入1次，总共输入745组数据。

图5 算例1地层岩质边坡坡脚处输入的地震加速度时程曲线

根据计算所得数据绘制各块体单元水平方向、竖直方向的加速度时程曲线。通过分析各个块体单元加速度的极值发现，块体单元1的水平方向、竖直方向加速度峰值与块体单元5的水平向、竖直向加速度峰值差值最大，说明地震波对这两个块体单元产生的影响最为显著。块体单元1和5的水平、竖直方向加速度时程曲线如图6、7所示。各块体单元水平方向、竖直方向加速度极值见表2。

分析图6、7可知，在地震载荷作用下，块体单元沿水平方向、竖直方向的加速度在t=0～1.68 s间的波动频率较低，较为稳定；在t=1.68～9.82 s间的波动频率较高，边坡在地震作用下受到较大冲击；在t=9.82 s之后，水平方向、竖直方向加速度波动频率开始逐渐变缓，趋于稳定。

图6 算例1岩质边坡块体单元1、5水平方向和竖直方向地震加速度时程曲线

由表2可看出，在水平方向上，块体单元1的加速度最大值为2.124 m/s2、最小值为-2.057 m/s2；块体单元5的加速度最大值为1.446 m/s2、最小值为-1.041 m/s2，两块体单元的水平方向加速度最大值差值达到0.678 m/s2，最小值差值达到1.016 m/s2。在竖直方向上，块体单元1的加速度最大值为0.708 m/s2、最小值为-0.686 m/s2；块体单元5的加速度最大值为0.482 m/s2、最小值为-0.467m/s2；最大值差值达到0.226 m/s2，最小值差值达到0.219 m/s2。显而易见，不同块体因为受到不同的内力与地震荷载作用，输出的水平方向和竖直方向的加速度时程曲线也不尽相同。滑动面与结构面差异越大，则数值差值越大。

表2 算例1各块体单元加速度极值 m/s2

依据线性数学规划模型求解算例1岩质边坡的稳定性动态安全系数，地震作用下的0～14.9 s时间段内，安全系数时程曲线间隔0.02 s输出一次，共得到745组数据。图7为计算得出的边坡稳定性安全系数时程曲线。

图7 算例1岩质边坡的稳定性安全系数时程曲线

分析图7可知，在地震开始作用后的t=0～1.6 s时间段内，地震作用对边坡稳定性的影响较小，此时间段内边坡安全系数波动规律与自然状态下静态安全系数波动规律相符合。在t=1.68～9.82 s时间段内，地震作用变强，边坡稳定性安全系数的波动频率与幅度显著增大，在t=3.84 s时安全系数到达最大值2.010；在t=4.08 s时安全系数到达最小值1.017。当t=9.82 s时，安全系数波动幅度开始逐渐放缓。从边坡安全系数的时程曲线可以清晰直观地看出地震作用下边坡的动力稳定性安全系数随时间变化的情况，通过对比发现边坡的安全系数时程曲线波动规律与加速度波动规律相符合。在本次地震作用下，边坡的安全系数时程曲线数值始终大于1，可以判断本算例边坡符合稳定性要求。

采用7度地震烈度时岩质边坡动力稳定性安全系数最小值已经非常接近1，为了更好地体现地震作用对于岩质边坡动力稳定性的影响，采用多地震烈度的水平向设计地震加速度代表值进行计算。8、9度地震烈度条件下边坡强度储备系数-频数直方图如图8所示。岩质边坡不同地震烈度条件下的安全系数时程曲线如图9所示，安全系数极值见表3。

图8 地震烈度为8、9度条件下算例1岩质边坡强度储备系数-频数直方图

图9 算例1岩质边坡不同地震烈度条件下安全系数时程曲线

表3 算例1不同地震烈度边坡稳定性安全系数极值表

由图9和表3可知：(1)地震作用会显著降低岩质边坡的稳定性和安全性，在数值上表现为安全系数最小值大幅度减小，与现有研究结果相同。(2)在8度地震烈度条件下745组安全系数数据中，安全系数小于1的有23组；9度地震烈度条件下745组岩质边坡安全系数小于1的有60组。(3)地震烈度越大，则安全系数时程曲线的波动幅度越大，岩质边坡的失稳可能性越高，说明高地震烈度条件下边坡的稳定性分析至关重要。

本文方法可以在较高地震烈度下准确计算出地震作用期间岩质边坡各个时间点的安全系数，并绘制出岩质边坡安全系数的时程曲线。为了比较，本文分别计算了在不同地震烈度条件下时程分析法求得的下限解与拟静力法所求得的安全系数，结果见表4。由表4可看出，时程分析法(本文方法)所求得的最小安全系数小于拟静力法所求得的安全系数，随着地震烈度的增大，二者之间的差距逐渐增大。这是因为本文方法考虑了地震特性和土体的动力特性，计算结果相较于拟静力法精度更高，能够较好地评价边坡的失稳情况。

表4 算例1不同地震烈度下两种方法求得最小安全系数表

本文方法不仅可以计算并绘制岩质边坡的安全系数时程曲线，还可以根据计算结果进行内力分析，绘制各滑面上的内力时程曲线，再根据内力时程曲线计算绘制各滑面的屈服函数时程曲线，从而确定岩质边坡的主要破坏模式。从算例1中选取了如图10所示的块体单元结构作为算例2进行内力计算与分析，地震加速度峰值按8度地震选取，取值为0.3g。滑面BC上的内力时程曲线如图11所示，内力峰值见表5。

图10 算例2块体单元结构面示意图

分析图11与表5可知，结构面BC上法向力最大值出现在3.84 s，最大值为5 845 N，方向为负；最小值出现在4.08 s，最小值为2 377 N，方向为负，显然，法向力时程曲线波动规律与所求得的安全系数时程曲线一致。切向力最大值出现在4.08 s，最大值为4 094 N，方向为负；最小值出现在3.84 s，最小值为1 982 N，方向为负，显然，切向力时程曲线波动规律与所求得的安全系数时程曲线相反。

表5 算例2块体单元滑面BC上的内力极值 N

图11 算例2边坡结构面BC上的内力时程曲线

通过计算出各个时步滑面上的内力，可以计算出各个时步滑面的屈服函数并绘制结构面屈服函数时程曲线。屈服函数的计算公式为：

f(j)=|Vk|-tanφk/KjNk-cklk/Kj

(16)

式中：Vk为结构面上的切向力，N；Nk为结构面上的法向力，N；lk为结构面k的长度，m；ck为对应结构面k上的凝聚力，N/m；φk为对应结构面k上的内摩擦角，(°)。

当f(j)=0时滑面屈服；f(j)

为了确定本算例的主要破坏模式，选取图10所示的滑面1(AB-BC-CD)和滑面2(HF-FC-CD)两种滑动趋势进行对比分析。滑面1和滑面2的屈服函数时程曲线如图12所示。

图12 算例2边坡滑面1、2屈服函数时程曲线(地震裂度8度)

在8度地震烈度条件下，从输入地震波到结束的14.9 s内每个结构面均得到了745组屈服函数的计算结果。由图12(a)可知，滑面1结构面CD的屈服函数值绝大部分时间都等于0,结构面基本上已经完全屈服；结构面BC和结构面CD的屈服函数值只有少部分小于0，该两个结构面大部分时间都已经屈服。由图12(b)可知，结构面HF与结构面FC的屈服函数值绝大部分时间仍然小于0，滑面还未屈服。显然，地震作用下结构面HF和结构面FC仍然具有较高的抵抗变形的能力，因此滑面2的滑动趋势并非本算例的主要破坏模式，而结构面AB、BC、CD基本上均已丧失了抵抗变形的能力，因此滑面1(AB-BC-CD)的破坏模式应为本算例的主要破坏模式。

(1)本文基于动力学理论、极限平衡理论和线性数学规划理论提出一种地震作用下岩质边坡的动力稳定性的时程分析新方法。

(2)本文方法使用块体单元离散岩体，采用拟动力法理论计算块体单元的地震时程加速度，基于极限平衡理论建立岩质边坡动力稳定性的线性数学规划模型；使用循环迭代法计算岩质边坡动力稳定性安全系数，获得边坡动力稳定性安全系数的时程曲线并对滑体结构面进行内力分析，找出其主要的破坏模式。

(3)本文方法理论严谨、概念明确、计算精度和效率高、工程应用简便，可将其应用于岩质边坡的动力稳定性计算领域，可为岩质边坡的抗震设计提供参考。

(4)本文在进行研究分析时，所研究的边坡其滑动面是已知的，这是一种比较理想的情况，在今后的研究中将运用本文方法对实际地震作用下滑面未知的边坡进行稳定性分析与研究。

水资源与水工程学报2021年2期